這段課程的主題是「如何在地圖上呈現資料」,包含以下概念:
- 相關的 R 套件
- 實際案例練習
- 地點標示
- 事件密度
- 空間內差
本課程使用的套件:
- ggplot
- ggmap
- gstat
- sp
November 18, 2016
require(ggmap) map <- get_map(location = 'Taiwan', zoom = 7) ggmap(map)
地圖的位置是透過 location 參數來指定,而 zoom 則是控制地圖的大小
tpe = c(lon=121.5197,lat=25.0356) map <- get_map(location=tpe, zoom=11) ggmap(map)
map <- get_map(location=tpe, zoom=11, language="zh-TW") ggmap(map)
map <- get_map(location=tpe, zoom=11, language="zh-TW", maptype="roadmap") ggmap(map)
map <- get_map(location=tpe, zoom=11, language="zh-TW", maptype="satellite") ggmap(map)
map <- get_map(location=tpe, zoom=11, language="zh-TW", maptype="toner-lite") ggmap(map)
中央氣象局 2011-2016 年發生的地震資料
earthquake <- read.csv("cwb_earthquake.csv", stringsAsFactors = F, encoding="UTF-8")
str(earthquake)
## 'data.frame': 3879 obs. of 7 variables: ## $ No : chr "\u00a4p\u00b0\u03f0\xec" "105099 " "\u00a4p\u00b0\u03f0\xec" "\u00a4p\u00b0\u03f0\xec" ... ## $ Time : chr "2016/11/1 \u00a4U\u00a4\xc8 12:41:00" "2016/10/29 \u00a4U\u00a4\xc8 07:49:00" "2016/10/29 \u00a4U\u00a4\xc8 02:44:00" "2016/10/29 \u00a4U\u00a4\xc8 02:08:00" ... ## $ Lon : num 122 122 122 122 121 ... ## $ Lat : num 23.8 24.2 24.7 24.4 23.5 ... ## $ depth : num 37 30.1 29.6 18.6 15.9 19.1 19.4 31.6 9.4 34.1 ... ## $ scale : num 3.9 4.3 3.6 3.6 3.6 3.3 5.2 4.2 3.7 4 ... ## $ location: chr "\u00aa\u1f6c\u00bf\u00a4\u00acF\u00a9\u00b2\u00abn\u00a4\xe8 25.1 \u00a4\u00bd\u00a8\u00bd (\u00a6\xec\u00a9\U000ea87d\u00ac\"| __truncated__ "\u00aa\u1f6c\u00bf\u00a4\u00acF\u00a9\u00b2\u00a5_\u00a4\xe8 21.4 \u00a4\u00bd\u00a8\u00bd (\u00a6\xec\u00a9\U000ea87d\u00ac\"| __truncated__ "\u00a9y\xc4\U0017f92cF\u00a9\u00b2\u00a6\xe8\u00a4\xe8 20.2 \u00a4\u00bd\u00a8\u00bd (\u00a6\xec\u00a9\U000f7ce5_\u00a5\u00ab"| __truncated__ "\u00a9y\xc4\U0017f92cF\u00a9\u00b2\u00abn\u00a4\xe8 38.2 \u00a4\u00bd\u00a8\u00bd (\u00a6\xec\u00a9\U000e9e44\U0017f92bn\u00b"| __truncated__ ...
require(dplyr)
# 新增「年份」欄位
earthquake$year <- substr(earthquake$Time,1,4)
# 挑出 2014~2016 的資料
eq1416 <- filter(earthquake, year %in% c("2014","2015","2016")) %>%
select(Time, Lon, Lat, depth, scale, year) # 省略不必要的欄位
str(eq1416)
## 'data.frame': 1754 obs. of 6 variables: ## $ Time : chr "2016/11/1 \u00a4U\u00a4\xc8 12:41:00" "2016/10/29 \u00a4U\u00a4\xc8 07:49:00" "2016/10/29 \u00a4U\u00a4\xc8 02:44:00" "2016/10/29 \u00a4U\u00a4\xc8 02:08:00" ... ## $ Lon : num 122 122 122 122 121 ... ## $ Lat : num 23.8 24.2 24.7 24.4 23.5 ... ## $ depth: num 37 30.1 29.6 18.6 15.9 19.1 19.4 31.6 9.4 34.1 ... ## $ scale: num 3.9 4.3 3.6 3.6 3.6 3.3 5.2 4.2 3.7 4 ... ## $ year : chr "2016" "2016" "2016" "2016" ...
map <- get_map(location="Taiwan", zoom=7, language="zh-TW") ggmap(map) + geom_point(aes(x=Lon,y=Lat,colour=scale), data=earthquake, alpha=0.3)
ggmap(map) + geom_point(aes(x=Lon,y=Lat,colour=scale), data=eq1416, alpha=0.3) + facet_grid(~year)
ggmap(map, extent="panel", maprange=F) + geom_density2d(data=eq1416, aes(x=Lon,y=Lat))
# 黑白地圖較適合疊上色塊
map <- get_map(location="Taiwan", zoom=7, language="zh-TW", color="bw")
# 如同前一張圖的地圖和等高線
ggmap(map, extent="panel", maprange=F) +
# 統計轉換:換成
stat_density2d(data=eq1416, aes(x=Lon,y=Lat, fill=..level.., alpha=..level..),
size = 0.1, bins = 100, geom = 'polygon') +
# 色階及圖說設定
scale_fill_gradient("Earthquake Density", low = "green", high = "red") +
scale_alpha(range = c(0.10, 0.30), guide = FALSE) +
# 等高線移到較上面的圖層
geom_density2d(data=eq1416, aes(x=Lon,y=Lat, alpha=0.6)) +
# 全圖主題設定
theme(axis.title = element_blank(), text=element_text(size = 12))
facet_wrap(~year)
有時候我們只有「點」的觀測,但是卻需要整個「面」的資料,例如:氣象觀測都是點狀的,但是數值模式需要的是「網格」資料,這時候就需要空間內插。
需要的套件:
gstatspmaptoolsrequire(gstat) require(sp) require(maptools)
# 設定資料座標
coordinates(earthquake) = ~Lon+Lat
# 設定網格範圍
x.range <- as.numeric(c(117.5, 124.5))
y.range <- as.numeric(c(20.5, 27.2))
# 製作網格點
grd <- expand.grid(Lon = seq(from = x.range[1], to = x.range[2], by = 0.5),
Lat = seq(from = y.range[1], to = y.range[2], by = 0.5))
# 設定網格資料座標
coordinates(grd) <- ~Lon+Lat
gridded(grd) <- TRUE
plot(grd, cex = 1.5, col = "grey") points(earthquake, pch = 1, col = "red", cex = 0.3)
網格準備完成,接著我們把地震規模(scale)內插到網格點上。
# 呼叫 IDW 函數進行空間內插
idw <- idw(formula=scale~1, locations=earthquake, newdata=grd, debug.level=0)
# 轉換資料格式以利繪圖
idw.output = as.data.frame(idw, stringsAsFactors=F)
names(idw.output)[1:3] <- c("lon", "lat", "value")
套件 stat 也包含 krige 函數,使用方法與 idw 類似,可以參考函數說明。
現在,我們可以用 geom_tile() 把網格資料疊在地圖上了。
ggmap(map, extent="panel", maprange=F) +
geom_tile(data=idw.output, aes(x=lon, y=lat, fill=value), alpha=0.5) +
scale_fill_gradient(low = "green", high = "red") +
theme(axis.title = element_blank(), text = element_text(size = 12))
網格資料的美觀程度,取決於網格的密度,但是網格越密,內插要算越久。以下是把網格間距從 0.5 度調降到 0.1 度的結果。
我們示範了三種在地圖上顯示資訊的方式:
使用的套件: